home *** CD-ROM | disk | FTP | other *** search
/ Aminet 6 / Aminet 6 - June 1995.iso / Aminet / gfx / 3d / irit50src.lha / irit5 / symb_lib / distance.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-09-23  |  18.8 KB  |  483 lines
  1. /******************************************************************************
  2. * Distance.c - Distance computations.                          *
  3. *******************************************************************************
  4. * Written by Gershon Elber, May 93.                          *
  5. ******************************************************************************/
  6.  
  7. #include <ctype.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include "symb_loc.h"
  11.  
  12. #define SELF_INTER_REL_EPS 100
  13.  
  14. static CagdPtStruct
  15.     *GlblSrfDistPtsList = NULL;
  16.  
  17. static void SymbSrfDistFindPointsAux(CagdSrfStruct *Srf,
  18.                      CagdRType Epsilon,
  19.                      CagdBType SelfInter);
  20. static void SrfDistAddZeroPoint(CagdRType u, CagdRType v, CagdRType Epsilon);
  21.  
  22. /*****************************************************************************
  23. * DESCRIPTION:                                                               M
  24. * Given a curve and a point, finds the nearest point (if MinDist) or the     M
  25. * farest location (if MinDist FALSE) from the curve to the given point.      M
  26. *   Returned is the parameter value of the curve.                 M
  27. *   Computes the zero set of (Crv(t) - Pt) . Crv'(t).                 M
  28. *                                                                            *
  29. * PARAMETERS:                                                                M
  30. *   Crv:        The curve to find its nearest (farest) point to Pt.          M
  31. *   Pt:         The point to find the nearest (farest) point on Crv to it.   M
  32. *   MinDist:    If TRUE nearest points is needed, if FALSE farest.           M
  33. *   Epsilon:    Accuracy of computation.                                     M
  34. *                                                                            *
  35. * RETURN VALUE:                                                              M
  36. *   CagdRType:   Parameter value in the parameter space of Crv of the        M
  37. *                nearest (farest) point to point Pt.                         M
  38. *                                                                            *
  39. * KEYWORDS:                                                                  M
  40. *   SymbDistCrvPoint, curve point distance                                   M
  41. *****************************************************************************/
  42. CagdRType SymbDistCrvPoint(CagdCrvStruct *Crv,
  43.                CagdPType Pt,
  44.                CagdBType MinDist,
  45.                CagdRType Epsilon)
  46. {
  47.     CagdRType TMin, TMax, t,
  48.     ExtremeDist = MinDist ? INFINITY : -INFINITY;
  49.     CagdPtStruct *TPt,
  50.     *Pts = SymbLclDistCrvPoint(Crv, Pt, Epsilon);
  51.  
  52.     /* Add the two end point of the domain. */
  53.     CagdCrvDomain(Crv, &TMin, &TMax);
  54.     TPt = CagdPtNew();
  55.     TPt -> Pt[0] = TMin;
  56.     TPt -> Pnext = Pts;
  57.     Pts = TPt;
  58.     TPt = CagdPtNew();
  59.     TPt -> Pt[0] = TMax;
  60.     TPt -> Pnext = Pts;
  61.     Pts = TPt;
  62.  
  63.     for (TPt = Pts, t = TMin; TPt != NULL; TPt = TPt -> Pnext) {
  64.     int i;
  65.     CagdPType EPt;
  66.     CagdRType
  67.         Dist = 0.0,
  68.         *R = CagdCrvEval(Crv, TPt -> Pt[0]);
  69.  
  70.     CagdCoerceToE3(EPt, &R, - 1, Crv -> PType);
  71.  
  72.     for (i = 0; i < 3; i++)
  73.         Dist += SQR(EPt[i] - Pt[i]);
  74.  
  75.     if (MinDist) {
  76.         if (Dist < ExtremeDist) {
  77.         t = TPt -> Pt[0];
  78.         ExtremeDist = Dist;
  79.         }
  80.     }
  81.     else {
  82.         if (Dist > ExtremeDist) {
  83.         t = TPt -> Pt[0];
  84.         ExtremeDist = Dist;
  85.         }
  86.     }
  87.     }
  88.     CagdPtFreeList(Pts);
  89.  
  90.     return t;
  91. }
  92.  
  93. /*****************************************************************************
  94. * DESCRIPTION:                                                               M
  95. * Given a curve and a point, find the local extremum distance points on the  M
  96. * curve to the given point.                                                  M
  97. *   Returned is a list of parameter value with local extremum.             M
  98. *   Computes the zero set of (Crv(t) - Pt) . Crv'(t).                 M
  99. *                                                                            *
  100. * PARAMETERS:                                                                M
  101. *   Crv:        The curve to find its extreme distance locations to Pt.      M
  102. *   Pt:         The point to find the extreme distance locations from Crv.   M
  103. *   Epsilon:    Accuracy of computation.                                     M
  104. *                                                                            *
  105. * RETURN VALUE:                                                              M
  106. *   CagdPtStruct *:   A list of parameter values of extreme distance         M
  107. *              locations.                                             M
  108. *                                                                            *
  109. * KEYWORDS:                                                                  M
  110. *   SymbLclDistCrvPoint, curve point distance                                M
  111. *****************************************************************************/
  112. CagdPtStruct *SymbLclDistCrvPoint(CagdCrvStruct *Crv,
  113.                   CagdPType Pt,
  114.                   CagdRType Epsilon)
  115. {
  116.     int i;
  117.     CagdCrvStruct *ZCrv,
  118.     *DCrv = CagdCrvDerive(Crv);
  119.     CagdPType MinusPt;
  120.     CagdPtStruct *ZeroSet;
  121.  
  122.     for (i = 0; i < 3; i++)
  123.     MinusPt[i] = -Pt[i];
  124.  
  125.     Crv = CagdCrvCopy(Crv);
  126.     CagdCrvTransform(Crv, MinusPt, 1.0);
  127.  
  128.     ZCrv = SymbCrvDotProd(Crv, DCrv);
  129.     CagdCrvFree(Crv);
  130.     CagdCrvFree(DCrv);
  131.  
  132.     ZeroSet = SymbCrvZeroSet(ZCrv, 1,  Epsilon);
  133.     CagdCrvFree(ZCrv);
  134.  
  135.     return ZeroSet;
  136. }
  137.  
  138. /*****************************************************************************
  139. * DESCRIPTION:                                                               M
  140. * Given a curve and a line, finds the nearest point (if MinDist) or the      M
  141. * farest location (if MinDist FALSE) from the curve to the given line.         M
  142. *   Returned is the parameter value of the curve.                 M
  143. *   Let Crv be (x(t), y(t)). By substituting x(t) and y(t) into the line     M
  144. * equation, we derive the distance function.                     M
  145. *   Its zero set, combined with the zero set of its derivative provide the   M
  146. * needed extreme distances.                             M
  147. *                                                                            *
  148. * PARAMETERS:                                                                M
  149. *   Crv:        The curve to find its nearest (farest) point to Line.        M
  150. *   Line:       The line to find the nearest (farest) point on Crv to it.    M
  151. *   MinDist:    If TRUE nearest points is needed, if FALSE farest.           M
  152. *   Epsilon:    Accuracy of computation.                                     M
  153. *                                                                            *
  154. * RETURN VALUE:                                                              M
  155. *   CagdRType:   Parameter value in the parameter space of Crv of the        M
  156. *                nearest (farest) point to line Line.                        M
  157. *                                                                            *
  158. * KEYWORDS:                                                                  M
  159. *   SymbDistCrvLine, curve line distance                                     M
  160. *****************************************************************************/
  161. CagdRType SymbDistCrvLine(CagdCrvStruct *Crv,
  162.               CagdLType Line,
  163.               CagdBType MinDist,
  164.               CagdRType Epsilon)
  165. {
  166.     CagdRType TMin, TMax, t,
  167.     ExtremeDist = MinDist ? INFINITY : -INFINITY;
  168.     CagdPtStruct *TPt,
  169.     *Pts = SymbLclDistCrvLine(Crv, Line, Epsilon);
  170.  
  171.     /* Add the two end point of the domain. */
  172.     CagdCrvDomain(Crv, &TMin, &TMax);
  173.     TPt = CagdPtNew();
  174.     TPt -> Pt[0] = TMin;
  175.     TPt -> Pnext = Pts;
  176.     Pts = TPt;
  177.     TPt = CagdPtNew();
  178.     TPt -> Pt[0] = TMax;
  179.     TPt -> Pnext = Pts;
  180.     Pts = TPt;
  181.  
  182.     for (TPt = Pts, t = TMin; TPt != NULL; TPt = TPt -> Pnext) {
  183.     CagdPType EPt;
  184.     CagdRType
  185.         Dist = 0.0,
  186.         *R = CagdCrvEval(Crv, TPt -> Pt[0]);
  187.  
  188.     CagdCoerceToE2(EPt, &R, - 1, Crv -> PType);
  189.  
  190.     Dist = EPt[0] * Line[0] + EPt[1] * Line[1] + Line[2];
  191.     Dist = ABS(Dist);
  192.  
  193.     if (MinDist) {
  194.         if (Dist < ExtremeDist) {
  195.         t = TPt -> Pt[0];
  196.         ExtremeDist = Dist;
  197.         }
  198.     }
  199.     else {
  200.         if (Dist > ExtremeDist) {
  201.         t = TPt -> Pt[0];
  202.         ExtremeDist = Dist;
  203.         }
  204.     }
  205.     }
  206.     CagdPtFreeList(Pts);
  207.  
  208.     return t;
  209. }
  210.  
  211. /*****************************************************************************
  212. * DESCRIPTION:                                                               M
  213. * Given a curve and a line, finds the local extreme distance points on the   M
  214. * curve to the given line.                             M
  215. *   Returned is a list of parameter value with local extreme distances.      M
  216. *   Let Crv be (x(t), y(t)). By substituting x(t) and y(t) into the line     M
  217. * equation, we derive the distance function.                     M
  218. *   Its zero set, combined with the zero set of its derivative provide the   M
  219. * needed extreme distances.                             M
  220. *                                                                            *
  221. * PARAMETERS:                                                                M
  222. *   Crv:        The curve to find its nearest (farest) point to Line.        M
  223. *   Line:       The line to find the nearest (farest) point on Crv to it.    M
  224. *   Epsilon:    Accuracy of computation.                                     M
  225. *                                                                            *
  226. * RETURN VALUE:                                                              M
  227. *   CagdPtStruct *:   A list of parameter values of extreme distance         M
  228. *              locations.                                             M
  229. *                                                                            *
  230. * KEYWORDS:                                                                  M
  231. *   SymbLclDistCrvLine, curve line distance                                  M
  232. *****************************************************************************/
  233. CagdPtStruct *SymbLclDistCrvLine(CagdCrvStruct *Crv,
  234.                  CagdLType Line,
  235.                  CagdRType Epsilon)
  236. {
  237.     CagdPType Pt;
  238.     CagdCrvStruct *TCrv1, *CrvX, *CrvY, *CrvZ, *CrvW, *DistCrv, *DerivDistCrv;
  239.     CagdPtStruct *ZeroSet1, *ZeroSet2, *TPt;
  240.  
  241.     SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ);
  242.     if (CrvZ)
  243.     CagdCrvFree(CrvZ);
  244.  
  245.     Pt[0] = Pt[1] = Pt[2] = 0.0;
  246.     CagdCrvTransform(CrvX, Pt, Line[0]);
  247.     CagdCrvTransform(CrvY, Pt, Line[1]);
  248.     TCrv1 = SymbCrvAdd(CrvX, CrvY);
  249.     CagdCrvFree(CrvX);
  250.     CagdCrvFree(CrvY);
  251.     if (CrvW) {
  252.     CagdCrvTransform(CrvW, Pt, Line[2]);
  253.     DistCrv = SymbCrvAdd(TCrv1, CrvW);
  254.     CagdCrvFree(CrvW);
  255.     CagdCrvFree(TCrv1);
  256.     }
  257.     else {
  258.     Pt[0] = Line[2];
  259.     CagdCrvTransform(TCrv1, Pt, 1.0);
  260.     DistCrv = TCrv1;
  261.     }
  262.  
  263.     ZeroSet1 = SymbCrvZeroSet(DistCrv, 1,  Epsilon);
  264.  
  265.     DerivDistCrv = CagdCrvDerive(DistCrv);
  266.     CagdCrvFree(DistCrv);
  267.  
  268.     ZeroSet2 = SymbCrvZeroSet(DerivDistCrv, 1,  Epsilon);
  269.     CagdCrvFree(DerivDistCrv);
  270.  
  271.     if (ZeroSet1 == NULL)
  272.     return ZeroSet2;
  273.     else if (ZeroSet2 == NULL)
  274.     return ZeroSet1;
  275.     else {
  276.     for (TPt = ZeroSet1; TPt -> Pnext != NULL; TPt = TPt -> Pnext);
  277.  
  278.     TPt -> Pnext = ZeroSet2;
  279.  
  280.     return ZeroSet1;
  281.     }
  282. }
  283.  
  284. /*****************************************************************************
  285. * DESCRIPTION:                                                               M
  286. * Given two curves, creates a bivariate scalar surface representing the      M
  287. * distance function square, between them.                     M
  288. *                                                                            *
  289. * PARAMETERS:                                                                M
  290. *   Crv1, Crv2:  The two curves, Crv1(u) and Crv2(v), to form their distance M
  291. *                function square between them as a bivariate function.       M
  292. *                                                                            *
  293. * RETURN VALUE:                                                              M
  294. *   CagdSrfStruct *:   The distance function square d2(u, v) of the distance M
  295. *                      from Crv1(u) to Crv2(v).                              M
  296. *                                                                            *
  297. * KEYWORDS:                                                                  M
  298. *   SymbSrfDistCrvCrv, curve curve distance                                  M
  299. *****************************************************************************/
  300. CagdSrfStruct *SymbSrfDistCrvCrv(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2)
  301. {
  302.     CagdSrfStruct *TSrf, *DiffSrf, *DotProdSrf,
  303.     *Srf1 = CagdPromoteCrvToSrf(Crv1, CAGD_CONST_U_DIR),
  304.     *Srf2 = CagdPromoteCrvToSrf(Crv2, CAGD_CONST_V_DIR);
  305.  
  306.     if (CAGD_IS_BSPLINE_SRF(Srf1) || CAGD_IS_BSPLINE_SRF(Srf2)) {
  307.     CagdRType UMin1, UMax1, VMin1, VMax1, UMin2, UMax2, VMin2, VMax2;
  308.  
  309.     if (CAGD_IS_BEZIER_SRF(Srf1)) {
  310.         TSrf = CnvrtBezier2BsplineSrf(Srf1);
  311.         CagdSrfFree(Srf1);
  312.         Srf1 = TSrf;
  313.     }
  314.     if (CAGD_IS_BEZIER_SRF(Srf2)) {
  315.         TSrf = CnvrtBezier2BsplineSrf(Srf2);
  316.         CagdSrfFree(Srf2);
  317.         Srf2 = TSrf;
  318.     }
  319.  
  320.     CagdSrfDomain(Srf1, &UMin1, &UMax1, &VMin1, &VMax1);
  321.     CagdSrfDomain(Srf2, &UMin2, &UMax2, &VMin2, &VMax2);
  322.     BspKnotAffineTrans(Srf1 -> VKnotVector,
  323.                Srf1 -> VLength + Srf1 -> VOrder,
  324.                VMin2 - VMin1, (VMax2 - VMin2) / (VMax1 - VMin1));
  325.     BspKnotAffineTrans(Srf2 -> UKnotVector,
  326.                Srf2 -> ULength + Srf1 -> VOrder,
  327.                UMin1 - UMin2, (UMax1 - UMin1) / (UMax2 - UMin2));
  328.     }
  329.  
  330.     DiffSrf = SymbSrfSub(Srf1, Srf2);
  331.     DotProdSrf = SymbSrfDotProd(DiffSrf, DiffSrf);
  332.     
  333.     CagdSrfFree(Srf1);
  334.     CagdSrfFree(Srf2);
  335.     CagdSrfFree(DiffSrf);
  336.  
  337.     return DotProdSrf;
  338. }
  339.  
  340. /*****************************************************************************
  341. * DESCRIPTION:                                                               M
  342. * Given a scalar surface representing the distance function square between   M
  343. * two curves, finds the zero set of the distance surface, if any, and        M
  344. * returns it.                                     M
  345. *   The given surface is a non negative surface and zero set is its minima.  M
  346. *   The returned points will contail the two parameter values of the two     M
  347. * curves that intersect in the detected zero set points.             M
  348. *                                                                            *
  349. * PARAMETERS:                                                                M
  350. *   Srf:         A bivariate function that represent the distance square     M
  351. *                function between two curves.                     M
  352. *   Epsilon:     Accuracy control.                                           M
  353. *   SelfInter:   Should be consider self intersection? That is, is Srf       M
  354. *                between a curve to itself!?                     M
  355. *                                                                            *
  356. * RETURN VALUE:                                                              M
  357. *   CagdPtStruct *:  A list of parameter values of both curves, at all       M
  358. *                    detected intersection locations.                        M
  359. *                                                                            *
  360. * KEYWORDS:                                                                  M
  361. *   SymbSrfDistFindPoints, curve curve distance, curve curve intersection    M
  362. *****************************************************************************/
  363. CagdPtStruct *SymbSrfDistFindPoints(CagdSrfStruct *Srf,
  364.                     CagdRType Epsilon,
  365.                     CagdBType SelfInter)
  366. {
  367.     GlblSrfDistPtsList = NULL;
  368.  
  369.     if (Srf -> GType == CAGD_SBEZIER_TYPE) {
  370.     Srf = CnvrtBezier2BsplineSrf(Srf); /* To keep track on the domains. */
  371.     SymbSrfDistFindPointsAux(Srf, Epsilon, SelfInter);
  372.     CagdSrfFree(Srf);
  373.     }
  374.     else
  375.     SymbSrfDistFindPointsAux(Srf, Epsilon, SelfInter);
  376.  
  377.     return GlblSrfDistPtsList;
  378. }
  379.  
  380. /*****************************************************************************
  381. * DESCRIPTION:                                                               *
  382. * Auxiliary function of SymbSrfDistFindPoints - does the hard work.         *
  383. *                                                                            *
  384. * PARAMETERS:                                                                *
  385. *   Srf:         A bivariate function that represent the distance square     *
  386. *                function between two curves.                     *
  387. *   Epsilon:     Accuracy control.                                           *
  388. *   SelfInter:   Should be consider self intersection? That is, is Srf       *
  389. *                between a curve to itself!?                     *
  390. *                                                                            *
  391. * RETURN VALUE:                                                              *
  392. *   None                                                                 *
  393. *****************************************************************************/
  394. static void SymbSrfDistFindPointsAux(CagdSrfStruct *Srf,
  395.                      CagdRType Epsilon,
  396.                      CagdBType SelfInter)
  397. {
  398.     int i, j,
  399.     ULength = Srf -> ULength,
  400.     VLength = Srf -> VLength;
  401.     CagdRType
  402.     *XPts = Srf -> Points[1];
  403.  
  404.     for (i = 0; i < ULength; i++) {
  405.     for (j = 0; j < VLength; j++) {
  406.         if (*XPts++ <= 0.0) {
  407.         CagdRType UMin, UMax, VMin, VMax;
  408.  
  409.         CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
  410.  
  411.         if (SelfInter) {
  412.             CagdRType
  413.             SelfInterEps = Epsilon * SELF_INTER_REL_EPS;
  414.  
  415.             if (UMax - UMin < SelfInterEps &&
  416.             VMax - VMin < SelfInterEps &&
  417.             UMax - VMin < SelfInterEps &&
  418.             VMax - UMin < SelfInterEps)
  419.                 return;
  420.         }
  421.  
  422.         /* The surface may have a zero here. Test the domain size   */
  423.         /* to make sure it makes sense to subdivide.            */
  424.         if (UMax - UMin < Epsilon && VMax - VMin < Epsilon) {
  425.             SrfDistAddZeroPoint((UMin + UMax) / 2.0,
  426.                     (VMin + VMax) / 2.0, Epsilon);
  427.         }
  428.         else {
  429.             /* Subdivide the surface and invoke recursively. */
  430.             CagdSrfStruct *Srf1, *Srf2;
  431.             CagdRType t;
  432.             CagdSrfDirType Dir;
  433.  
  434.             if (UMax - UMin > VMax - VMin) {
  435.             t = (UMin + UMax) / 2.0;
  436.             Dir = CAGD_CONST_U_DIR;
  437.             }
  438.             else {
  439.             t = (VMin + VMax) / 2.0;
  440.             Dir = CAGD_CONST_V_DIR;
  441.             }
  442.             Srf1 = CagdSrfSubdivAtParam(Srf, t, Dir);
  443.             Srf2 = Srf1 -> Pnext;
  444.             SymbSrfDistFindPointsAux(Srf1, Epsilon, SelfInter);
  445.             SymbSrfDistFindPointsAux(Srf2, Epsilon, SelfInter);
  446.             CagdSrfFree(Srf1);
  447.             CagdSrfFree(Srf2);
  448.         }
  449.         return;
  450.         }
  451.     }
  452.     }
  453. }
  454.  
  455. /*****************************************************************************
  456. * DESCRIPTION:                                                               *
  457. * Auxiliary function of SymbSrfDistFindPoints - part of the hard work.         *
  458. *                                                                            *
  459. * PARAMETERS:                                                                *
  460. *   u, v:   A zero set location found to be added to global parameter list.  *
  461. *   Epsilon:    To match against existing zero set point, for simiarity,     *
  462. *                                                                            *
  463. * RETURN VALUE:                                                              *
  464. *   void                                                                     *
  465. *****************************************************************************/
  466. static void SrfDistAddZeroPoint(CagdRType u, CagdRType v, CagdRType Epsilon)
  467. {
  468.     CagdPtStruct *Pt;
  469.  
  470.     for (Pt = GlblSrfDistPtsList; Pt != NULL; Pt = Pt -> Pnext) {
  471.     if (ABS(Pt -> Pt[0] - u) < Epsilon * 10 &&
  472.         ABS(Pt -> Pt[1] - v) < Epsilon * 10)
  473.         return;                     /* Point is already there. */
  474.     }
  475.  
  476.     Pt = CagdPtNew();
  477.  
  478.     Pt -> Pt[0] = u;
  479.     Pt -> Pt[1] = v;
  480.     Pt -> Pnext = GlblSrfDistPtsList;
  481.     GlblSrfDistPtsList = Pt;
  482. }
  483.